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An  indirect  (pointer-based)  datastructure  for  th^/ storage  of  a  symmetric,  positive  definite 
matrix  of  coefficients,  and  the  corresponding  “right  hand  side”  vector  is  presented.  This 
datastructure  is  useful  in  many  kinds  of  engineering  work,  most  notably  finite  element  analysis. 
Procedures  for  accessing  the  datastructure  are  demonstrated  via  the  small  Turbo  Pascal 
program  KxR.  * 


2.  Introduction 


Matrix  methods  in  general,  and  the  finite  element  method  (FEM)  in  particular,  are  now 
common  engineering  tools,  y  The  equations  to  be  solved  in  the  vast  majority  of  finite  element 
stress  analyses  are  — — _ _ _  _ _  * 
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or,  in  matrix  notation, 


Kx 


where  K  is  the  positive  definite}  stiffness  matrix,  x  is  the  vector  of  unknowns,  and  R  is  the 
vector  of  knowns.  In  the  general  case,  of  course,  the  number  of  equations  is  much  larger  than 
6,  often  numbering  into  the  thousands. 


The  majority  of  techniques  for  solving  (1)  are  variants  of  a  method  introduced  by  the 
mathematical  genius  Carl  F.  Gauss.  See,  for  example,  [Ralston]}  or  [Carnahan],  for  an 
explication  of  the  technique.  The  demonstrations  in  [Bathe]  and  [Wilson]  are  outstanding. 
The  original  Gaussian  method  has  been  supplanted  in  recent  times  by  methods  that  are 
appropriate  for  computer  implementation. 

2.1  Motivation 

Since  there  are  published  programs,  why  bother  writing  yet  another  solverT  The  major  need 
was  for  a  solver  of  “reasonable”  capacity,  in  Pascal,  for  a  personal  computer  (PC). 

t.1.1  Pascal  m  general  Pascal,  a  programming  language  invented  in  the  early  1970’s  by  Niklaus 
Wirth,  is  intrinsicly  legible  code.  That  is,  a  reader  and  writer  of  the  code  can  understand  a 
Pascal  program  more  easily  than  code  written  in  other  languages.  Pascal  also  facilitates  data 
structuring,  a  language  capability  that  is  absent  from  FORTRAN  and  BASIC.  Finally,  Pascal  is  a 
readable  step  along  the  pathway  to  C,  the  inscrutable  (but  ubiquitous)  language  in  which  it 
seems  most  programs  will  eventually  be  written. 

t.l.t  Turbo  Pateal  m  particular  Version  5.0  of  the  Turbo  Pascal  compiler  manufactured  by 
Borland  International  [BorlandPU] ,  [BoriandR]  provides  a  superb  built-in  debugger.  Borland 
has  implemented  the  Modula-2  module  concept,  which  allows  separate  compilation  of  code 
modules  while  retaining  the  strong  type-checking  of  Pascal.}  The  module  concept,  present  in 


t  A  matrix  Q  is  pet***  itfinUt  if  the  product  xrQx  >  0  for  any  non-iero  X. 
t  References  we  listed  st  the  end  of  the  paper, 
f  Borland’s  spelling  for  the  word  “module”  is  “Unit” 
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almost  all  languages  except  BASIC  and  standard  Pascal,  speeds  up  program  development 
enormously.  Overlays  are  available  in  Turbo  Pascal  to  allow  very  large  programs  to  fit  into  the 
DOS3x  memory  limitation  of  640K  bytes.  And,  finally,  Turbo  Pascal  is  faster  than  any  other 
compiler  for  the  PC,  for  any  language,  especially  FORTRAN.  This  extremely  high  speed  makes 
for  a  good  environment  for  program  development,  for  “small”  programs. 

£.1.8  Drawback »  to  Turbo  Pascal  Turbo  Pascal  is  net  Pascal:  some  of  Turbo  Pascal  is  not 
portable  to  other  machines. 
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S.  Date  Storage  on  PC* 


3.1  Capacity 

Hie  maximum  addressable  storage  on  a  DOS  machine  is  usually  abbreviated  to  the  number 
640K  bytes.  The  architecture  of  the  80x86  series  machines  accesses  storage  in  tegmenta.  The 
maximum  size  of  a  single  segment  is  65520  bytes. 

3.2  Floating  Point  Numbers 

Turbo  Pascal  allows  five  representations  of  floating  point  numbers,  as  shown  in  Table  1: 


|  Table  1:  Floating  Point 

in  Turbo  Pascal 

TVpe 

4 

16380 

real 

6 

10920 

double 

8 

8190 

extended 

10 

6552 

comp 

10 

6552 

The  first  column  gives  the  name  of  the  Pascal  floating  point  type,  and  the  second  column  shows 
the  number  of  bytes  required  for  that  type.  The  third  column  of  the  table  shows  the  maximum 
number  of  floating  point  variables  allowed  in  a  segment  This  effectively  limits  the  maximum 
size  of  data  a  program  can  use,  and  the  maximum  size  of  a  single  array,  unless  other  techniques 
are  used.  Some  PC  FORTRAN  compilers  map  arrays  to  appropriate  portions  of  segments,  but 
control  of  this  process  is  not  available  in  the  language  direedy. 

Interestingly  enough,  computations  are  faster  using  the  Turbo  Pascal  “double”  than  the 
“real”,  (using  a  numerical  co-processor  chip)  because  the  truncation  of  the  result  of  floating 
point  arithmetic  operations  is  avoided.  Certainly,  computations  are  more  accurate  in  “double” 
than  in  “real”  or  “single”. 


4.  A  Comparison  of  Matrix  Storage  Techniques 

The  crucial  datastructure  in  a  FEM  program  is  the  stiffness  matrix.  Because  the  majority  of  the 
procedures  not  involved  with  reading  input  or  printing  output  are  related  to  accessing  this  great 
datastructure,  the  form  of  the  stiffness  matrix  dictates  the  structure  of  the  FEM  code. 

Usually,  one  wishes  to  mesh  a  structure  with  the  economically  largest  possible  number  of 
elements  ( numel)  because  accuracy  of  analysis  is  directly  related  to  the  fineness  of  the  grid. 
Imagine  a  mesh  of  2D  planar  8-noded  isoparametric  elements.  These  elements  have  two 
degrees  of  freedom  per  node,  and  therefore,  16  equations  per  element.  In  order  to  solve  a  2D 
FEM  mesh  with  ONE  of  these  elements,  one  must  be  able  to  store  and  solve  at  least  16 
equations.  Let  us  consider  the  following  four  storage  forms  for  K  Square,  Triangular,  Banded, 
and  Skyline. 


4.1  Square 

An  n-dimensional  square  matrix  may  be  illustrated  with 
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where  n=»8  in  this  case.  Assuming  that  calculations  will  be  performed  using  “doubles”,  the 
total  storage  for  this  matrix  is  n2  <  8190  where  n  is  the  sise  of  K  Thus,  the  maximum  n  is 
90,  and  the  upper  bound  on  numel  in  a  2D  mesh  is  about  8,  which  is  interesting,  but  not  very 
useful. 

A  complete  (and  naive)  Pascal  program  to  solve  equation  (l)  is  developed  in  [Wirth]  using  the 
stepwise  refinement  technique. 
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program  gauss  .elimination  ; 
coast 

n  =6  ; 
var 

i,  j,  k:  1  ..  n  ; 

p  :  real ; 

A  :  array  [  1  ..  n,  1  ..  n  ]  of  real ; 

B  :  array  [  1  ..  n  ]  of  real ; 

X  :  array  [  1  ..  n  ]  of  real  ; 

begin 

{  assign  values  to  A  and  B  } 
for  k  :=1  to  n  do  begin 
p  :=1.0  /  A[k,k  ]  ; 
for  j  :=k  +  1  to  n  do 

A[k,j  ]  p  *  A [k, j  ]  ; 

B[k  ]  :=p  *  B[k  ]  ; 
for  i  :=k  +  1  to  n  do  begin 
for  j  :=k  +  1  to  n  do 

A(i,j  ]  :=A[i,j  ]  -  A[i,k  ]  *  A[kfj  ]  ; 

B(i  ]  :=*B[i  ]  -  A[i,k  ]  *  B[k  ] 
end  { i  } 
end  {  k  }  ; 
k  :=n  ; 
repeat 

p:-B[k]  ; 

for  j  :*=k  +  1  to  n  do 

p  :=*p  -  A[k,j  ]  *  X[ j  ]  ; 

X(k  ]  :=p  ; 
k  k  -  1 
until  k  =0 

{X[l  ]  ...  X[n  ]  are  the  solutions  } 

end. 

Program  Is  Gaussian  Solution  of  Simultaneous  Equations 

Program  1  gives  the  flavor  of  Gaussian  elimination,  but  needs  some  improvement.  There 
should  at  least  be  a  check  for  division  by  sero. 

4.2  Triangular 

Because  K  is  symmetric,  or  can  be  made  so  by  appropriate  multiplications,  only  "half”  the 
array,  either  above  or  below  the  diagonal,  need  be  stored. 
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The  storage  required  for  K  is  n(»+l)/2  <  8190  ,  so  n  <  127  ,  and  the  maximum  nttmei  is 
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about  12.  His  is  an  improvement  over  the  square  scheme,  but  hardly  enough  to  make  the 
derivation  of  a  structure  to  store  the  array  worthwhile. 

4.3  Banded 


A  further  observation  of  FEM  stiffness  equations  shows  that  for  judiciously  numbered  meshes, 
the  matrix  appears  to  be  “banded”;  ie.,  all  of  the  non-zero  terms  in  K  are  within  a  certain 
distance  from  the  diagonal.  E.g., 
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If  it  is  known  that  certain  values  of  K  are  zero,  then  they  need  not  be  stored.  This  leads  to  the 
auxiliary  descriptor  of  K  the  maximum  bandwidth,  commonly  represented  in  programs  as  the 
variable  mazbw.  The  smaller  the  bandwidth,  the  less  storage  required. 


array  of  the  form 
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where  the  asterisks  indicate  “wasted”  space.  Storage  required:  n-maxbw  <  8190  .  Thus,  n 
varies  from  90  (with  a  mazbw  of  90)  to  8190  (with  a  mazbw  of  1).  Recalling  the  simplistic  2D 
mesh,  the  maximum  number  of  elements  is  8190/18  —=511.  Now  we’re  getting  somewhere! 
Of  course,  the  algorithm  that  solves  equation  (l)  when  Kia  stored  in  banded  form  will  be  more 
complicated  than  Program  1,  because  it  will  have  to  map  the  subscripts  of  K  to  the  coordinates 
of  the  rectangular  array. 

4.4  Skyline 

Wilson,  and  others,  have  observed  that  K  stored  in  banded  form  may  store  unnecessary  zeroes. 
The  bandwidth  at  any  location  in  K  may  vary  from  small  to  large,  forming  a  “skyline”,  as 
shown  typically  below  for  n  —17.  Note  that  the  solvers  described  in  [Wilson]  and  [Bathe] 
eliminate  computations  that  use  terms  outside  the  skyline,  and  are  among  the  fastest  solvers  in 
common  use. 
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Since  the  “columns”  are  of  arbitrary  height,  so  a  little  more  information  about  K  than  n  and 
maxbw  is  needed.  Unfortunately,  FORTRAN  has  no  direct  facility  for  storing  this  datastructure. 


The  technique  used  in  [Wilson]  and  [Bathe]  to  store  the  skyline  array  in  the  programs  SAP4, 
NONSAP  and  ADINA  conceptually  numbers  each  element  of  the  array  consecutively.  The 
indexing  scheme  proceeds  from  the  diagonal  element,  numbered  first,  through  each  element  in 
the  column,  numbering  towards  the  skyline,  viz., 
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K  is  then  stored  in  a  one-dimensional  (linear)  array,  with  an  auxiliary  integer  array  D  to  point 
to  the  diagonal  terms  of  K,  thus: 


i  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18 

D[i]  1  2  4  6  9  11  13  15  17  26  29  33  37  41  44  55  58  71 

For  example,  the  indexes  of  all  terms  in  the  upper  portion  of  column  9  are  defined  by  the 
limits  D[9]  to  D[10]- 1  (—  17  to  25).  The  last  term  (D[18])  is  asentinel  pointer  to  the  n+lth 
column.  Wilson  uses  these  FORTRAN  “pointers”  to  calculate  column  heights,  and  requires  an 
additional  pointer  to  enable  calculation  of  the  height  of  the  last  column. 
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The  profile  P  is  defined  mathematically  as 

P  =  EDf 

iM 

which  is  simply  the  number  of  non-zero  terms  within  the  skyline.  The  total  storage  (in  bytes) 
required  to  store  K  must  now  include  storage  for  D: 

Storage  =  PSizeof(float)  +  ( n +1) -Sizeof(  integer) 

where  Sizeof  returns  the  number  of  bytes  required  to  store  its  argument.  The  skyline  technique 
will  probably  be  a  gain  over  the  banded  storage  scheme  although  K  is  still  stored  in  a  single 
segment. 

Note  that  the  trade-offs  made  to  gain  storage  space  will  make  die  solution  algorithm  more 
complex  and  harder  to  understand  than  Program  1.  Since  Wilson  uses  FORTRAN,  that 
language’s  built-in  features  to  clearly  calculate  subscripts  must  be  supplemented  by  array  access 
mechanisms  which  map  K  to  the  skyline.  Now  consider  Pascal,  a  language  that  has  the  facility 
for  direct  access  to  unusual  and  arbitrary  datastructures. 
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S.  Implementation  of  the  Skyline  Data  Structure  in  T\irbo  Pascal 

So  far,  about  the  best  that  any  of  the  storage  schemes  yield  is  an  upper  limit  on  the  number  of 
2D  elements  somewhere  around  500,  certainly  enough  for  programs  analyzing  small  models.  Is 
it  possible  to  do  any  better  without  resorting  to  a  file-based  storage  scheme?  The  answer  is  of 
course  “yes”,  but  it  requires  that  the  use  of  pointers,  a  data  structuring  tool  present  in  most 
modern  languages. 

5.1  Pointers 

Prom  [Jensen],  Chapter  10: 

“A  static  variable  (staticly  allocated)  is  one  which  is  declared  in  a 
program  and  subsequently  denoted  by  its  identifier.  It  is  called  static,  for 
it  exists  (  ie.,  memory  is  allocated  for  it)  during  the  entire  execution  of 
the  block  to  which  it  is  local.  A  variable  may,  on  the  other  hand,  be 
generated  dynamically  (without  any  correlation  to  the  static  structure  of 
the  program)  by  the  procedure  n ew.  Such  a  variable  is  consequently 
called  a  dynamic  variable. 

Dynamic  variables  do  not  occur  in  an  explicit  variable  declaration  and 
cannot  be  referenced  directly  by  identifiers.  Instead,  generation  of  a 
dynamic  variable  introduces  a  pointer  value  (which  is  nothing  other  than 
the  storage  address  of  the  newly  allocated  variable).  Hence,  a  pointer 
type  P  consists  of  an  unbounded  (tic)  set  of  of  values  pointing  to 
elements  of  a  given  type  T.  P  is  then  said  to  be  bound  to  T.  The  value 
nil  is  always  an  element  of  P  and  points  to  no  element  at  all. 


type  <identifier>  =  *  <type  identifier  > 


If,  for  example,  p  is  a  pointer  to  a  variable  of  type  T  by  the  declaration 

var  p  :  *  T 

then  p  is  a  reference  to  a  variable  of  type  T,  and  p"  denotes  that 
variable.!  In  order  to  create  or  allocate  such  a  variable,  the  standard 
procedure  new  is  used.  The  call  new(p)  allocates  a  variable  of  type  T  and 
assigns  its  address  to  p. 

Pointers  are  a  simple  tool  for  the  construction  of  complicated  and  flexible 
data  structures.  If  the  type  T  is  a  record  structure  that  contains  one  or 
more  fields  of  type  ‘T,  then  structures  equivalent  to  arbitrary  finite 
graphs  may  be  built,  where  the  Ps  represent  the  nodes,  and  the  pointers 
are  the  edges.” 

Pointers  have  several  advantages:  they  permit  programmers  to  develop  the  most  flexible  of  data 
structure;  they  permit  fast  access  to  data  items  without  need  for  subscript  computation;  they  do 
not  require  recourse  to  assembly  language;  and  they  allow  access  to  all  of  available  memory. 
The  latter  is  crucial.  For  example,  the  memory  available  to  a  PC  program  is  typically  470 
Kbytes,  (approximately  58750  doubles  as  opposed  to  the  8190  indicated  in  Table  1).  By  using 


t  If  p  is  a  pointer  to  a  type  X  the  expression  p*  is  said  to  dereference  the  pointer.  Alternatively,  p  is  the  address,  p*  is 
the  "valne”  at  the  address. 
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pointers  to  access  this  storage,  the  upper  bound  on  n time/  for  our  simple  2D  problem  becomes 
5874  elements!  Furthermore,  the  upper  bound  on  name/  for  a  20-node  3D  isoparametric 
element  mesh  is  1631.  The  pointer  data  storage  technique  therefore  affords  the  capability  to 
solve  significant  FEM  problems  in  main  storage  on  the  PC. 

The  First  Law  of  Thermodynamics  paraphrased  states  that  you  don’t  get  something  for  nothing. 
Great  flexibility  means  that  the  programmer  must  impose  discipline  to  control  code  and  data. 
(Pointers  have  been  called  the  datastructure  equivalent  of  the  GOTO  statement.)  The  pointer 
concept  is  foreign  to  the  FORTRAN  programmer,  so  a  significant  number  of  engineers  are 
unlikely  to  easily  understand  pointer  code  without  sufficient  time  to  study  the  concept.  The 
program  becomes  more  intellectually  dense,  since  more  concepts  are  embedded  in  each 
character.  Access  to  a  data  item  is  by  indirection,  rather  than  directly.  Both  programmer  and 
code  readers  not  only  move  their  lips  when  reading  pointer  code,  but  read  out  loud,  assume 
awkward  postures,  and  draw  little  diagrams.  So  much  for  “intrinsicly  legible  code”  in  Pascal. 
However,  pointers  are  supported  by  the  language  itself,  so  learning  and  using  pointers  becomes 
more  natural  and  straight-forward  with  practice. 

5.2  Memory  allocation  In  T\irbo  Pascal 

The  memory  available  to  a  Turbo  Pascal  programmer  comes  from  the  “Heap”,  an  area  in 
memory  accessible  by  several  Turbo  procedures  and  functions. 

Dispose  Releases  storage  previously  allocated  to  a  variable.  Used  in  cooperation  with 

New. 

Freemem  returns  to  the  heap  exactly  the  number  of  bytes  allocated  to  a  pointer  by 
Getmem. 

Getmem  allocates  any  number  of  bytes  (and  therefore,  partial  records)  to  a  pointer. 

Used  in  cooperation  with  Freemem. 

Mark  establishes  a  pointer  to  a  location  on  the  heap;  used  in  cooperation  with 

Release.  . 

Memavail  returns  the  memory  available  at  the  current  point  in  the  program’s  execution. 

Maxavail  returns  the  maximum  memory  available  on  the  heap. 

New  allocates  storage  to  a  variable.  Used  in  cooperation  with  Dispose. 

Release  frees  all  memory  above  a  pointer  established  by  Mark. 

The  cooperating  procedure  pairs  are  therefore:  New,  Dispose;  Mark,  Release;  and  Getmem, 
Freemem.  Getmem  and  Freemem  are  non-standard  Pascal,  and  are  the  closest  to  the  C 
language  storage  allocation  routines. 

5.3  The  New  Data  Structure 
The  ingredients  include: 

1.  Memavail  from  Turbo  Pascal.  This  procedure  will  indicate  the  total  storage  available. 

2.  Mark  and  Release  from  Turbo  Pascal.  Mark  will  allocate  the  total  heap  storage  to  the 
arrays  needed  for  the  solution  of  (1).  Release  will  return  the  heap  storage  to  the 
operating  system  on  completion  of  the  job.  Although  this  may  seem  unnecessarily 
fastidious,  my  machine  has  run  out  of  memory  during  the  development  of  KxR,  because 
“no  one  told  him”  that  he  regained  the  memory. 

3.  Getmem  from  Turbo  Pascal  (the  key  ingredient).  Getmem  will  be  used  to  allocate  to 
each  column  only  that  storage  needed  to  eontain  the  column. 
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4. 


height,  a  pointer  to  a  variable  length  integer  array  that  describes  the  height  of  each 
column.  This  replaces  the  D  array  of  Wilson’s  method,  and  all  the  explicit  calculations  of 
columnar  bounds  are  handled  directly  by  the  language. 

5.  K,  a  pointer  to  an  array  of  pointers  containing  the  elements  of  K.  The  stiffness  matrix  is 
therefore  accessed  doubly  indirectly  (pointers  to  pointers).  Although  this  may  seem 
complicated,  it  actually  makes  the  code  simpler  and  clearer  than  the  linear  FORTRAN 
technique. 


The  skyline  structure  in  Pascal  will  therefore  be  represented  in  a  manner  similar  to  Wilson’s. 
However,  each  column  of  the  array  will  be  numbered  from  the  diagonal  to  the  uppermost 
element,  rather  than  numbering  each  element  of  the  array.  This  alternative  can  be 
demonstrated: 


Figure  lx  Conceptual  Layout  of  K 


The  array  of  column  heights  is  then 

i  123456789  10  11 

height*[i]  12232222934 


12  13  14  15  16 
4  4  3  11  3 


where  we  write  “height*[ij”  to  show  what  the  va/ves  are.  Remember,  “height”  is  just  a 
pointer;  the  presence  of  the  '**”  dereferences  the  pointer. 
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With  this  groundwork  laid,  it  is  possible  to  write  the  code  to  implement  this  structure. 
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8.  The  code 


The  allocation  of  the  data  structure  is  discussed  first,  followed  by  the  solver  itself.  Additional 

utility  routines  for  access  to  the  data  structure  follow. 

8.1  Allocation  at  Storage 

The  allocation  of  storage  follows  this  algorithm: 

1.  Assume  the  number  of  equations  neq  is  known. 

2.  Then  an  array  of  neq  integers  is  needed  to  hold  height  Allocate  this  array. 

3.  Determine  the  height  of  each  column  and  store  it  in  the  appropriate  location  of  height. 

4.  Since  neq  and  height  are  known,  the  shape  of  K  is  completely  determined.  Allocate 
storage  for  an  array  of  neq  pointers  to  columns  of  K  (as  i  varies  from  1  to  neq,  K*[t] 
points  to  column  t).  To  each  pointer  K‘[t],  allocate  AetjAf  *[t]  floating  point  values.  With 
this  scheme,  the  term  on  the  diagonal  in  the  i-th  column  is  written  K*[t]  ‘[l]  and  the  j-th 
term  from  the  diagonal  in  the  t-th  column  is  written  K*[ »']  *[/]. 

5.  Allocate  neq  floating  point  values  for  the  right  hand  side  vector  R  . 

The  Turbo  Pascal  procedures  to  accomplish  the  algorithm  outlined  above  are  shown  below. 
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6.1.1  Allocation  of  height  Hie  data  types  involved  are 


TYPE 

float 

=  double  ;  ( *  change  to  suit  accuracy  *) 

CONST 

MAXEQNS 

=  65520  div  Sizeof(  float )  ; 

TYPE 

column  Jieight_array 

=  array  (  1  ..  MAXEQNS  ]  of  word  ; 

column  Jieightjtype 

=  ‘column Jieight_array  ; 

VAR 

height : 

columnjieightjype  ; 

Note  that  Turbo  Pascal  allows  the  declarative  matter  of  a  program  to  be  in  arbitrary  order  (not 
standard  Pascal).  Also,  simple  computations  are  permitted  in  the  declarations:  efiv  is  the  Pascal 
“integer  divide”  operation.  The  type  word  permits  variables  of  that  type  to  assume  values 
between  0  and  65532.  Float  is  defined  as  double  only  once  in  the  solver;  thus  the  accuracy  may 
be  changed  with  the  alteration  of  this  single  switch.  MAXEQN  is  arbitrarily  defined  as  the 
maximum  number  of  floats  that  will  fit  in  a  single  memory  segment. 

The  procedure  to  allocate  height: 

( ***** 

* 

*.  allocate  Jieight  -  allocate  storage  for  the  column  heights. 

* 

*) 

procedure  allocate_height ; 

VAR 

c  :  longint ; 

begin  (*  allocate  Jieight  *) 

e  :==  neq  *  sizeof(  word  )  ; 

if  c  >  MemAvail  then  begin 

writeln(  ’Out  of  memory  attempting  to  allocate  ’,  c, 

’  bytes  to  column  heights  -  allocate  Jieight.’  )  ; 

abort 

end 

else 

Getmem  (  height,  c  ) 
end  ( *  allocate  Jieight  *)  ; 

In  addition  to  the  procedures  already  mentioned,  allocate Jieight  makes  use  of  abort: 
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( ***** 

* 

*.  abort  -  release  memory  and  die. 

* 

*) 

procedure  abort ; 
begin  ( *  abort  *) 

release(  origin  )  ; 
halt(  1  ) 

end  ( *  abort  *)  ; 

which  presumes  that  the  initialization  of  the  program  includes  the  statement 
mark(  origin  )  ; 

Halt  is  the  Turbo  Pascal  procedure  which  terminates  program  execution. 
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6.1.2  Allocation  of  K  Declarations  of  the  data  types  associated  with  the  storage  of  K  and  R  : 


TYPE 

column  Jype  —  array  [  1  ..  MAXEQNS  ]  of  float ; 

column_ptr  =  *column_type  ; 


column_array  =  array  (  1  ..  MAXEQNS  ]  of  column _ptr  ; 
matrix  =  *column_array  ; 


VAR 

K  :  matrix  ; 

R  :  column_ptr  ; 


The  procedure  to  allocate  the  storage  for  K  begins  with  “allocate _JC”  which  allocates  storage  to 
the  pointers  to  the  columns: 

( ***** 

* 

*.  allocate _JC  -  allocate  storage  for  the  matrix  of  coefficients. 


procedure  allocate _JC  ; 

VAR 

c  :  longint ; 
i  :  word  ; 

begin  (*  allocate_JC  *) 

c  :=  neq  *  sizeof(  pointer  )  ; 

if  c  >  MemAvail  then  begin 
writeln  (  ’Out  of  memory  attempting  to  allocate  c, 
’  bytes  to  [K]  -  allocate JC.’  )  ; 

abort 

end 

else 

Getmem  (  K,  c  )  ; 

for  i  :=-=  1  to  neq  do 
allocatejcolumn(  i,  height*[i]  ) 

end  (*  allocate  JC  *)  ; 
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The  first  part  of  this  procedure  is  similar  to  allocate _Jieight.  In  the  second  part,  the  for 
statement  calls  allocatejcolumn  to  allocate  the  storage  for  each  column  of  K.  This  procedure 
is: 

^  ***** 

* 

*.  allocate_column  -  arrange  storage  for  the  i-th  column. 

* 

*) 

procedure  allocate_column  (  i :  integer  ;  ( *  the  column  of  interest  *) 
r  :  word  ( *  the  size  desired  *) 

)  ; 

VAR 

c  :  longint ; 
begin 

c  :=  r  *  sizeof_fioat ; 

if  c  >  MemAvail  then  begin 
writeln  (  ’Out  of  memory  attempting  to  allocate  c, 

’  bytes  to  column  i,  ’  -  allocatejcolumn.’  )  ; 

abort 

end 

else  begin 

Getmem  (  K*[i] ,  c  )  ; 
end 

end  ( *  allocatejcolumn  *)  ; 

Allocatejcolumn  uses  a  previously  computed  tizeof_jioat,  which  is 
siseof_fioat  :=*=  Si*  eof(  float); 

Note  that  the  neat  arrangement  of  K  in  Figure  1  is  not  necessarily  how  the  datastructure  is 
arranged  in  memory.  Getznem  searches  the  heap  for  the  next  block  of  the  appropriate  size,  and 
assigns  it  to  the  pointer.  Because  the  /angaage  accommodates  this  process,  the  actual  locations 
in  storage  assigned  to  each  column  of  K  are  of  little  concern. 


6.1.S  Allocation  of  R  Hie  right  hand  side  vector  is  allocated  by  the  following: 

[***** 

* 

*.  allocate_vector  -  allocate  storage  for  a  ’right-hand-side- like’  vector. 
*) 

procedure  allocate.vector  (  VAR  R  :  column_ptr  )  ; 

VAR 

c  :  Ion  gin  t ; 

begin  (*  allocate.vector  *) 

c  :=  neq  *  siieof Jloat ; 

if  c  >  MemAvail  then  begin 
writeln  (  ’Out  of  memory  attempting  to  allocate  ’,  c, 

’  bytes  to  column  vector  -  allocate_vector.’  )  ; 

abort 

end 

else 

Getmem  (  R,  c  ) 
end  (*  allocate_vector  *)  ; 

All  storage  necessary  to  the  solution  of  (1)  has  now  been  allocated. 
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0.2  The  solver 

The  solver  itself  is  an  adaptation  of  the  LDl  r  technique  demonstrated  in  [Wilson]  and  [Bathe] . 
Control  of  the  solution  is  via  the  executive  procedure  active_column_solver: 

( ***** 

* 

*.  active_columnj5olver  -  Gaussian  solution  of  simultaneous  equations. 

*) 

procedure  active_column_solver ; 

begin  ( *  active_column_solver  *) 

write(  ’Decomposition  ...* )  ; 
decompose  ; 
write  In  (  ’  done.’ ); 

write(  ’Reduction  ...’ )  ; 
reduce  _jhs  ; 
writeln(  ’  done.’  )  ; 

write(  ’Back  substitution  ...’ )  ; 
back_substitute  ; 
writeln(  ’  done.’  ) 

end  ( *  active_column_polver  *)  ; 

Program  2at  Active  Column  Solver  -  Executive 

Messages  surround  the  major  phases  of  the  solution  process  so  the  solution  may  be  monitored 
from  the  screen. 
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6.2.1  Decomposition  Hie  decomposition  is  handled  by  decompose: 

*.  decompose  •  LDU  decomposition  of  coefficient  matrix. 

* 

*) 

procedure  decompose  ; 

VAR 

b,  c  :  float ; 

n,  P>  <1>  *»  t,  u,  v  :  word  ; 
begin  ( *  decompose  *) 
for  n  1  to  neq  do  begin 

(* 

*  propagate  effect  of  previous 

*  off-diagonal  terms  to  the  current  column 
*) 

if  heighten]  >  1  then  begin 
(*  p  -=  first  previous  column  of  interest  *) 
p  n  -  heighten]  +  2  ; 
for  q  :=  heighten]  -  1  down  to  2  do  begin 
s:»2; 
t  :=  q  +  1  ; 

(*  u  —  number  of  multiplies  to  perform  *) 
u  min(  h eight*  [p]  -  1,  heighten]  —  q )  ; 
c  0.0  ; 

for  v  :==  1  to  u  do  begin 
c:-c+  K*[p]‘[s]  *  K*[n]  *[t]  ; 

Inc(  s  )  ; 

Inc(  t ) 
end  (*  v  *)  ; 

K*[n]‘[q]  :-K>]*(q]-  c; 

Inc(  p  )  (*  the  next  ”previous  column”  *) 
end  (*  q  *) 
end  (*  if  *)  ; 

(* 

*  propagate  effect  of  previous 

*  diagonal  terms  to  the  current  column 
*) 

p n  ; 

b  0.0  ; 

for  q  :*=  2  to  heighten]  do  begin 
Dec(  p  )  ; 

c  K‘[nj  *(q]  /  K*(p]  *(l]  ;  (*  the  ”Li f  terms  *) 

b:-b+  c  *K*[n]‘[q]  ; 

K*[n]*[q]:=c 
end  (*  q  *)  ; 

K*[n]‘(l]  :-K*[nj*[l]-  b; 
if  K‘[n]*[l]  <-*0.0  then  solerr(  n  ) 
end(*n*) 

end  (*  decompose  *)  ; 

Program  2bt  Active  Column  Solver  -  Decomposition 
Although  the  “*”  notation  may  look  a  little  strange,  the  code  bean  a  striking  resemblance  to 
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non-pointer  code  (compare  with  Program  1).  The  pointers  obviate  explicit  calculation  of 
subscripts  of  terms  in  K,  however,  so  the  code  is  relatively  clean.  While  studying  decompose, 
bear  in  mind  that  K*(p]‘[l]  is  the  diagonal  term  in  column  p,  and  K‘[p]*[2]  is  the  first  off- 
diagonal  term.  Probably,  the  most  efficient  way  to  observe  how  the  solver  works  is  to  prepare  a 
small  test  case  (similar  to  the  example  demonstrating  use  of  KxR),  and  watch  the  program  run 
with  the  Turbo  Debugger. 

The  function  min  is  simply: 

^  ***** 

* 

*.  min  -  returns  smallest  integer  value. 

* 

*) 

function  min  (  a,  b  :  integer  )  :  integer  ; 
begin  ( *  min  *) 

if  a  <=  b  then  min  a  else  min  :=  b 
end  ( *  min  *)  ; 


Inc  and  Dec  are  Turbo  Pascal  procedures  that,  respectively,  increment,  and  decrement,  their 
arguments.  They  map  directly  to  machine  instructions,  so  are  quite  fast. 

In  case  of  trouble,  decompose  calls  solan 
( ***** 

*.  sol  err  -  prints  error  message  and  dies. 

* 

*) 

procedure  solerr  (  n  :  integer  )  ; 
begin  (*  solerr  *) 

writeln(  ****  Diagonal  term  in  column  n,  ’  ■*=  ’,  K*[n]  “[l]  )  ; 
write  In  (  *  Coefficient  matrix  is  NOT  positive  definite.  -  solerr.’ )  ; 
halt(  1  ) 

end  (*solerr*); 

Solerr  should  probably  use  abort,  rather  than  halt. 
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S.S.t  Reduction  The  reduction  of  the  right  hand  aide  is  performed  by  reduce_yhs: 

(  ***** 

* 

*.  reducejhs  -  using  factors  computed  by  'decompose’. 

* 

*) 

procedure  reduce_jhs  ; 

VAR 
c  :  float ; 
n,  p,  q  :  word  ; 
begin  (*  reducejhs  *) 
for  n  :=  1  to  neq  do  begin 
p  :=  n  -  1  ;  (*  first  "previous  column"  *) 

c  :=  0.0  ; 

for  q  :=  2  to  heighten]  do  begin 
c  :=  c  +  K*[n]*[q}  *  R*[p]  ; 

Dec(  p  )  (*  next  "previous  column”  *) 

end  (  *  q  *)  ; 

R*[n]  :=  R*[n]  -  c 
end  ( *  n  *)  ; 
end  (*  reduce  j-hs  *)  ; 

Program  2a  Active  Column  Solver  -  Reduction 
6. 2.8  Back  Subetitution  The  back  substitution  step  is  handled  by  bsekjrabartitute: 

^  ***** 

* 

*.  backjiubstitute  •  complete  solution  of  the  equations. 

* 

*) 

procedure  backjiubstitute  ; 

VAR 

n,  p,  q  :  word  ; 
begin  (*  backjiubstitute  *) 
for  n:*l  to  neq  do 
R*[n]  :=*=  R*[n]  /K‘[n]*[l]  ; 
for  n  neq  downto  2  do  begin 
P  n  -  1  ; 

for  q  :*  2  to  heighten]  do  begin 
R*[p]  :*=  R*[P]  ■  K*[n]  *[q]  •  R*(n]  ; 

Dec(  p ) 
end  (*  q  *) 
end  (*  n  *) 

end  (*  backjubstitute  *)  ; 

Program  2d  Active  Column  Solver  -  Back  Substitution 
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0.3  Accessing  the  Data  Structure  vis  Utility  Routines 

Several  of  the  procedures  and  functions  that  provide  an  interface  to  the  datastructure  are 
commented  on  below.  A  more  complete  set  is  listed  in  the  Appendix. 

6.S.1  Clearing  K  In  order  to  use  the  datastructure,  most  FEM  programs  will  require  a  “clean 
slate.”  The  following  procedure  sets  the  contents  of  Kto  0.0: 

^ ***** 

* 

*.  clear _f£  -  set  the  contents  of  K  to  0.0. 

* 

*) 

procedure  clear JK  ; 

VAR 
i :  word  ; 

begin  ( *  clear  JK  *) 
for  i  :=  1  to  neq  do 
set_vector(  K*(i] ,  height*[i] ,  0.0  ) 
end  ( *  clearJC  *)  ; 

CJear_K  makes  use  of  the  more  general  procedure  set_vectcr  (below).  The  danger  lurking 
here  is  that  undisciplined  use  of  aet_yector  will  clear  more  storage  than  has  been  allocated, 
possible  destroying  the  program  or  operating  system  code. 

( ***** 

* 

*.  set_vector  -  set  V*[l..n]  to  W. 

*) 

procedure  setvector  (  V  :  column_ptr;  n  :  word;  W  :  float )  ; 
begin  (  *  set_vector  *) 
for  n  :=  1  ton  do 
V*[n]  :=  W 
end  ( *  set_vector  *)  ; 
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6.S.2  Adding  a  Term  to  K  Given  the  location  at_row,  at_col  in  a  “full”  K,  this  procedure  adds 
term  to  the  appropriate  location  in  K  stored  in  skyline  form.  That  is,  add_to_K  maps  the 
standard  matrix  subscripts  to  the  skyline.  Relatively  cheap  insurance  is  provided  by  the  if 
statements,  which  avoid  storing  a  term  in  a  random  memory  location. 

( ***** 

* 

*.  add_to_K  -  add  ’term’  at  ’at_row,at_col’ 

*  where  at_row,  at_col  give  the  location  in 

*  a  "full*  matrix. 

* 

*) 

procedure  add  to  K  (  term  :  float ;  at_row,  at_col  :  word  )  ; 

VAR 
i  :  word  ; 

begin  (*  add_to_K  *) 

if  at_row  <=  at_col  then  begin  (*  K  is  upper  triangle  only  *) 
i  :=  Succ(  at_col  -  atjrow  )  ; 
if  i  <=  height“[at_col]  then  (‘insurance*) 

K*[at_col]  *[i]  :=  K*[at_col]  *(i]  +  term 
end  ( *  if  *) 
end  ( *  addjto _JC  *)  ; 

6.S.S  Vie  Product  f  =*  K  *  x  The  usual  definition  of  a  matrix  multiplying  a  vector  is  f  =  Kx 

■ 

or,  using  indices,  /,•  =  Let  K  be  symmetric,  and  consider  the  equation 

/- » 


/ 1 

*11  *12  *13 

f  ' 

*1 

/  2 

ft 

.  . 

*12  *22  *23 
*13  *23  *33 

*2 

*3 

»  4 

Performing  the  multiplication  yields  the  following  terms  for  f: 

/ 1  *=*11*1  +  *18*2  +  *13*3 
/  a  =  *12*1  +  *22*2  +  *28*3 
/ 3  ***13*1  +  *23*2  +  *33*8 

The  datastructure  only  stores  “half”  of  K,  however.  Hence,  the  computation  of  the  product  is 
reordered  such  that  all  the  terms  on  and  above  the  diagonal  are  used  first,  thus: 

/ 1  =*ii*i 
/  2  =*12*1  +  *22*2 
ft  =*13*1  +  *83*2  +  *33*3 

This  of  coune  is  equivalent  to  using  all  the  terms  to  the  left  of  the  diagonal  in  a  full  matrix. 

Then,  each  column  is  considered  in  turn;  the  terms  at  the  appropriate  “row”  (=  height)  are 
used  to  accumulate  the  missing  terms  into  f:  Using  column  2, 

/ 1  =/ 1  +  *12*2 
and  using  column  3, 

/ 1  =/ 1  +  *13*3 
ft  =/ 2  +  *23*3 

It  is  seen  that  the  resulting  f  is  correct.  This  algorithm  is  implemented  by  the  following 
procedure.  Note  that  the  numbering  of  K  in  the  procedure  agrees  with  that  of  the 
datastructure. 
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( ***** 


*.  KtimesX  -  column  vector  F  :=  skyline  matrix  K  times  column  vector  X. 

*  uses  global  variables  ’neq’,  ’K’  and  ’height* 

* 

*  K  is  symmetric: 

*  kll  k44 

*  k21  k43 

*  k31  k42 

*  (sym)  k41 

* 

*  Part  1  uses  each  column  from  the  diagonal  to  the  maximum 

*  height  in  each  column  as  a  row.  This  ensures  that 

*  the  terms  to  the  left  of  the  diagonal  are  used  in  the 

*  multiplication. 

* 

*  Part  2  adjusts  the  previous  sums  by  using  the  terms  above 

*  the  diagonal  as  the  remainder  of  the  row. 

* 

*) 

procedure  KtimesX  (  F,  X  :  column_ptr  )  ; 

VAR 

i,  j.  q  :  word  ; 
s  :  float ; 

begin  (*  KtimesR  *) 
for  i  :=  1  to  neq  do  begin 
(*  Parti  •) 
q  height*  [i]  ; 

s  :=  0.0  ; 

for  j  :=  1  to  height*[i]  do  begin 
s  :=  s +■  K-[i]*[q]  *X*[i]  ; 
dec(  q  ) 
end  ; 

F*(i]  ==  s  ; 

(♦Part 2  *) 

q  :==  i  *  1  ; 

for  j  :=«  2  to  height*[i]  do  begin 
F*[q]  :==  F*[q]  +  K*[i]  *[j]  *  X*[i]  ; 
dec(  q  ) 
end 

end  ( *  for  *)  ; 
end  (*  KtimesX*); 


7.  Using  the  KxR  Program 

7.1  Command  line 

kxr  input  output 

7.2  Input  file 

The  input  is  free  format.  Each  entry  separated  by  blank,  tab,  or  <cr>.  The  maximum 
number  of  numbers  on  a  line  will  usually  depend  on  the  operating  system,  and  most  likely, 
requires  that  a  line  be  less  than  256  characters. 

It  is  good  practice  to  begin  each  kind  of  data  on  a  new  line. 

7.2.1  Input  Fit  Format  First,  on  a  line  by  itself,  the  number  of  equations. 

Second,  the  height  array.  This  may  use  as  many  lines  as  necessary  to  completely  describe  the 
array. 

Third,  the  terms  of  K,  in  order  from  the  diagonal  towards  the  top  of  each  column. 

Fourth,  the  terms  of  R,  from  first  equation  to  last. 

For  example,  consider  the  set  of  simultaneous  equations 


5.0  -4.0  1.0  0.0 

o 

-4.0  6.0  -4.0  1.0 

*3 

1 

1.0  -4.0  6.0  -4.0 

*3 

►  — —  ' 

0 

0.0  1.0  -4.0  5.0 

■  • 

*4 

0 

k  4 

Pertinent  data  for  this  set: 

—  The  number  of  equations  ntq  is  4. 

—  The  column  heights  are,  respectively,  1,  2,  3,  and  3. 

—  The  terms  of  column  1:  5.0. 

—  The  terms  of  column  2:  6.0,  and  -4.0. 

—  The  terms  of  column  3:  6.0,  -4.0,  and  1.0 

—  The  terms  of  column  4:  5.0,  -4.0,  and  1.0 

—  Finally,  the  terms  of  the  right  hand  side  are:  0.0,  1.0,  0.0,  and  0.0 

The  input  file  for  this  case  might  therefore  look  like: 

4 

12  3  3 


5.0 

6.0 

-4.0 

6.0 

-4.0 

1.0 

5.0 

-4.0 

1.0 

0.0 

1.0 

0.0 

7.3  Output  file 

Hie  output  file  contains  the  solution  vector,  one  term  per  line.  Since  the  solution  vector  will 
most  likely  be  used  by  another  program,  there  is  no  label  or  “pretty-print”  information  written 
to  this  file. 

To  use  the  solution  vector  R  in  another  program: 

1.  Allocate  storage  for  ntq  floats,  to  hold  R. 

2.  Open  the  file  access  to  the  text  file  containing  R.  Turbo  Pascal  statements  would  have  the 
form 

var 

f  :  text ; 

assign (  f,  ’solution.dat’  )  ; 
reset(  f  )  ; 

3.  Read  R  with  statements  of  the  form 

for  i  :=  1  to  neq  do 
read(  f,  R*[i]  )  ; 
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APPENDIX 

The  following  pages  contain  the  complete  source  code  for  KxR. 
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float . Inc 


type 

float  -  doubla  ;  (•  change  to  suit  accuracy  •) 


••••••••••  global*. pas  ****•*****••**••••*••••••***•••• 

Unit  globala  ; 

(•.  globala  -  where  all  the  global  variable!  are  kept.  •) 


Interface  (*  public  declarations  •) 

{$1  float. ino  }  (*  defines  the  numerical  precision  •) 

CONST 

HAXEQNS  -  69320  div  Sizeof(  float  )  ; 


TYPE 

coluan_height_array «  array  [  1  ..  HAXEQNS  )  of  word  ; 
coluan_height_type  -  *coluan_height_array  ; 


coluan_type  ■  array  {  1  . .  HAXEQNS  ]  of  float  ; 

coluan_ptr  »  *coluan_type  ; 


coluan_array  «  array  (  1  ..  HAXEQNS  ]  of  coluan _ptr  ; 

aatrix  m  'coluan_array  ; 


aatrix.fi le 


«  file  of  float  ; 


itea  m  record  case  boolean  of 

true  :  (  w  :  word  )  ; 
false:  (  f  :  float  ) 
end  (•  itea  •)  ; 


VAR 


ooluan  :  coluan_ptr  ;  (•  pointer  to  the  K  coluan  *) 

f  :  text  ;  (•  input  file  *) 

g  :  text  ;  (•  output  file  •} 


height  :  coluan_height_type  ; 


K  :  aatrix  ; 

neq  :  word  ; 

afile  :  aatrix.file  ; 
afile.naae  :  string[l2] 

origin  :  pointer  ; 


(•  aatrix  of  coefficients  •) 

(•  »  of  equations,  this  problea  *) 


(•  narks  the  start  of  storage  *) 


pre.daax,  pre.dain  :  float  ; 
post_daax,  post.dain  :  float  ; 


R  :  eoluan_ptr  ; 


(•  "Right  hand  side”  of  equation  •) 


•  izao f_f loat  :  lonfint  ; 

X  :  coluan_ptr  ;  (•  used  to  chock  accuracy  •  ) 


Iaploaontation  (•  private  daclarationa  •) 

begin  (*  global*  initialization  coda  •) 

height  :»  nil  ; 

K  :»  nil  ; 

r  :»  nil  ; 

end  (•  global*  *) 


kxr.pt 


{*>+  > 

prograa  KxR  ; 


KxR  -  solve  the  aatrix  equation  [K] {x}  ■  {R> 
where  [R]  is  symmetric  and  stored  in 
skyline  fora  per  Bathe/Wilson. 

Copyright  (c)  Peter  N  Roth  -  January  1989 . 


Uses  overlay 
,  global* 
,  kxrutil 
.  init 
,  oalo 
.  fini 


{SO  kxrutil  > 

{SO  init  } 

{SO  oalo  > 

{SO  fini  } 

begin 

ovrinit(  'kxr.ovr'  )  ; 
initialise  ;  (•  unit  IK1T  •) 

oaloulate  ;  (•  unit  CALC  •) 

finish  (•  unit  F mi  •) 

end  (•  KxR  •) . 
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init.pi 


Unit  init  ; 

(*.  init  -  set  up  the  each in*.  read  input,  ate.  •) 

Interface  (•  public  declarations  •) 

Uses  overlay 
,  global* 

,  kxrutil 

(90  kxrutil  } 
procedure  initialize  ; 

Iepleaentation  (•  private  declaration*  •) 

VAR 

■X,  *2  :  longint  ;(•  teaporary  variable*  for  storage  calculations  •) 


(*•••• 

* 

*.  usage  -  'how  to  invoke  the  prograa.' 

« 

•) 

procedure  usage  ; 
begin  (*  usage  •) 

writeln(  '•••  usage:  paraastr(  O  ),  '  input  output  '  ) 

end  (•  usage  •)  ; 


*.  coaaand_line  -  get  and  vet  coaaand  line  paraaeters. 

•) 

procedure  coaaand_line  ; 
begin  (•  coaaand.line  *) 

if  paraacount  <•  x  then  begin 

usage  ; 

halt (  1  ) 

end  ; 

if  paraastr(  1  )  -  parses tr(  2  )  then  begin 

usage  ; 

writeln(  '  Input  and  output  file  naaes  aust  differ.'  )  ; 
halt (  l  ) 

end  ; 

assign (  f,  paraastr(  1  )  )  ; 
reset (  f  )  ; 

if  eof(  f  )  then  begin 

usage  ; 
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writeln(  '  laaediate  and  of  fllo  on  ',  paraastr(  1  )  )  ; 

usage  ; 

abort 

and  ; 

assign (  g,  paraastr(  2  )  )  ;  (*  will  overwrite  an  existing  file  *) 

rewrite (  g  ) 

and  (*  coaaand_line  •)  ; 


«.  get.neq  -  read  the  number  of  aquations  froa  <f>. 

•) 

procedure  get_neq  ; 
begin  (•  get_neq  •) 

readln(  f.  neq  )  ; 

if  neq  >  NAXEQNS  then  begin 

writeln(  'Nuaber  of  equations  exeeeds  current  capacity  (' , 

NAXEQNS ,  ')  -  get_neq. '  )  ; 

halt  (  l  ) 
end 

else  if  neq  <■  0  then  begin 

writeln(  'Nuaber  of  equations  <■  0  ???  Can"t  solve.  -  get_neq.'  )  ; 

abort 
end  (*  if  •) 

end  (•  get.neq  •)  ; 


«.  get_heights  -  read  the  heights  of  each  coluan. 


procedure  get_heights  ; 

VAR 

i  :  word  ; 

begin  (•  get_heights  •) 

for  i  :■  1  to  neq  do  begin 
read(  f.  height* fi]  )  ; 
if  height* [i]  <«  0  then  begin 

writeln(  'Non-positive  coluan  height  at  entry  ',i  )  ; 
abort 

end  (•  if  *) 
end  (•  for  •) 

end  (•  get .heights  •)  ; 
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«.  get_K  -  raad  the  coefficient  aatrix  K. 

•> 

procadura  get_K  ; 

VAR 

i,  i  :  word  ; 

bafin  (•  get_K  •) 

for  i  :•  1  to  naq  do  bagin 

for  j  :■  1  to  height* [i]  do  begin 
read<  f.  R*ti]*[j]  ) 

and 

and  (•  for  •) 
and  (•  gat_K  *)  ; 


*.  get_R  -  read  the  right  hand  aide  R. 
•) 

procadura  get_R  ; 

VAR 

i  :  word  ; 
bagin  (•  gat_R  •) 

for  i  l  to  naq  do 
read(f ,  R*[i]  ) 

and  (•  get_R  •)  ; 
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initialise  -  par fora  the  initialisations. 


• 
a. 
a 
•) 

procedure  initialise  ; 
begin  (•  initialise  •) 

coaaand_line  ;  (•  get  k  vet  the  coaaand  line  *) 

■1  : -  aeaavail  ; 

sizeof_float  :■  sizeof(  float  )  ; 

aark(  origin  )  ;  (•  where  the  heap  begins  *) 

get_neq  ; 

allocate_height  ; 
get_heights  ; 

(• 

•  allocate  ALL  space  first,  to  ensure 

•  that  we  can  coapute,  THEN  read 

•  the  stuff. 

•) 

allocate_K  ; 
allocate_vector(  R  )  ; 
allooate_vector(  X  )  ; 

(• 

•  It  is  not  necessary  to  clear  the  arrays. 

•  since  we  will  re-initialize  the  arrays 

•  via  the  READ. 

•) 

clearjt  ; 

set_vector(  R,  neq,  0.0  )  ; 
set_vector(  X.  neq,  0.0  )  ; 

gat_K  ; 
get_R  ; 

■2  :■  aeaavail  ; 

write (  'Meaory  available  (bytes):  al  )  ; 

writeln(  used:  ',  al  -  a2  )  ; 

assign (  afile,  'junk. eat'  )  ; 
rewrite (  afile  )  ; 

end  (•  initialise  *)  ; 

(• - + - ♦ - ♦ - ♦ - ♦ - - - •) 

end  (*  init  •) . 
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Unit  calc  ; 


(*.  calc  -  tha  aaat  of  tha  crunching.  •) 

Interface  (*  public  declaration*  *) 

Uses  overlay 
,  global* 

.  kxrutil 
,  solver 


(SO  kxrutil  > 

{$0  solver) 

procedure  calculate  ; 

Implementation  (•  private  declarations  *) 


(••••« 

•.  calculate  -  control  the  crunching. 

•) 

prooedure  calculate  ; 
begin  (•  calculate  *) 

diagonal_extreaes (  pre_daax,  pre_dain  )  ; 

active_ooluan_solver  ;  (•  unit  SOLVER  •) 

diagonal_extreaes (  post_daax,  post_dain  )  ; 

end  (*  calculate  •)  ; 


end  (•  calc  *) . 


f  ini.  pi 


Unit  fini  ; 

(*.  fini  -  clean  up  ft  go  hone.  •) 

Interface  (*  public  declarations  *) 

Uses  overlay 
,  globals 
,  kxrutil 


{$0  kxrutil  } 
procedure  finish; 

Implementation  (•  private  declarations  •) 


*.  write_R  -  put  the  solution  vector  to  <g>. 

•) 

procedure  write_R  ; 

VAR 

n  :  word  ; 

begin  (•  write_R  •) 

for  n  :■  1  to  neq  do 

writeln(  g,  R*[nJ  )  ; 

end  (•  write_R  *)  ; 


•.  finish  -  out  of  here. 

•) 

procedure  finish  ; 
begin  (•  finish  *) 

write_R  ; 

writeln(  'Diagonal  extremes  before  decomposition:'  )  ; 
writeln(  pre_dmax:30,  pre_dmin: 30  )  ; 

writeln(  'Diagonal  extremes  after  decomposition: '  )  ; 
writeln(  post_dmax:30,  post_dmin:30  }  ; 

release (  origin  )  ;  (•  memory  baok  to  system  •) 

close (  f  )  ; 
close(  g  ) 

end  (*  finish  *)  ;(•  — —  ♦ - +  — - —  +  — —  + - + - - - •) 

end  (•  fini  •). 


••••*•••••  solver. pas 


{*>♦} 

Unit  solver  ; 

(*.  solvsr.  •) 

Interface  (*  public  declarations  •) 

Uses  overlay 
.  globals 
,  kxrutil 


procedure  active_coluan_solver  ; 
procedure  decoapose  ; 
procedure  reduce_rhs  ; 
procedure  back_substitute  ; 


Implementation  (•  private  declarations  •) 
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decompose  -  LOU  decomposition  of  coefficient  matrix. 


•) 

procedure  decoapoae  ; 

VAX 

b,  c  :  float  ; 

n,  p,  q,  a,  t,  u,  v  :  word  ; 


basin  (•  decoapoae  •) 

for  n  :■  1  to  naq  do  basin 

(• 

•  propagate  affect  of  previous 

•  off-diagonal  teraa  to  the  currant  coluam 

•) 

if  height* [n]  >  1  than  begin 

(•pa  first  pravioua  column  of  interest  *) 

p  :«  n  -  height" [n]  +  3  : 

for  q  :■  height* Jn]  -  1  downto  2  do  begin 
a  :■  2: 
t  q  ♦  1  : 

(•  u  *  number  of  aultiplies  to  perform  •) 
u  ain(  height* [p]  -  1,  height* (n)  -  q  )  ; 
c  :■  0.0  ; 

for  v  :•  1  to  u  do  begin 

c  c  ♦  K*!pl*[s]  •  K*[n]*[t]  ; 

Inc(  a  )  ; 

Inc(  t  ) 
end  (•  v  •)  ; 

K*(n]*[q]  :•  K*[n]*[ql  -  c  ; 

Inc(  p  )  (•  the  neat  "previous  column"  •) 

end  (•  q  •) 
end  (•  if  •)  ; 

(* 

•  propagate  effect  of  previous 

•  diagonal  terms  to  the  current  column 

•) 

p  :•  n  ; 
b  0.0  ; 

for  q  :>  2  to  height* [n]  do  begin 
Dee(  p  )  ; 

c  K*tn]*t<l]  /  K*[P]*111  I  (•  the  "LiJ"  terms  •) 
b  :■  b  +  c  *  K*{nj*(qj  ; 

K*[n]*lq]  c 
end  (•  q  •)  ; 

K*ln]*[l]  K*tn]*[l]  -  b  : 

if  K*tn]*[l)  <•  0.0  then  solerr(  n  ) 
end  (•  n  •) 


end 


(•  deooapose  •)  ; 


*.  reduce_rhs  -  using  factors  coaputed  by  'dacoaposa' . 

*) 

procedure  reduce_rhs  ; 

VAR 

c  :  float  ; 
n  :  word  ; 
p  :  word  ; 
q  :  word  ; 

begin  (•  reduce_rhs  •) 

for  n  :■  1  to  neq  do  begin 

p  :■  n  -  1  ;  (•  first  "previous  coluan"  • 

e  :■  0.0  ; 

for  q  :■  2  to  height* [n)  do  begin 
o  c  ♦  K*[nriq]  •  R*[p]  ; 

Dec (  p  )  (•  next  "previous  coluan"  •) 

end  (•  q  •)  ; 

R*[nJ  a* In]  -  c 
end  («  n  •)  ; 

end  (•  reduce_rha  •)  ; 


*.  back_substitute  -  coaplete  solution  Of  the  equations. 
•) 

procedure  baek_aubstitute  ; 

VAR 

n  :  word  ; 
p  :  word  ; 
q  :  word  ; 

begin  (•  back_substitute  •) 

for  n  :■  1  to  neq  do 

R*[n]  R*tn]  /  R*[n]*ll]  I 

for  n  :«  neq  downto  2  do  begin 
P  n  -  1  ; 

for  q  :■  2  to  height* [n]  do  begin 
R'lPl  :•  **IP]  -  R‘[n]*[ql  •  R’[n)  5 
Deo(  p  ) 
end  (•  q  •) 
end  (•  n  •) 

end  (•  baok_substitute  •)  ; 
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*.  active_coluan_solvsr  -  Gaussian  solution  of  siaultanaous  aquations . 
•) 

procadura  set ivs_colu»n_solvsr  ; 

bagin  (•  activa_coluan_solvar  •) 

writs (  'Daeoaposition  ..."  )  ; 

daooaposa  ; 
writslnf  '  dona."  ); 

writs (  'Xaduction  ..."  )  ; 

raduoa_rhs  ; 
writsln(  '  dona."  )  ; 

writs (  'Back  substitution  ..."  )  ; 
back_substituta  ; 
writsln(  '  dona."  ) 

and  (•  act ira_coluan_solvar  •)  ; 
and  (•  solvsr  •) . 


ssssssssss  kxrutil.pas  ssssssssssssssssssssssssssssss 

{«o+} 

Unit  kxrutil  ; 

(*.  kxrutil  -  utility  routinas  for  tha  KxX  paokaga.  •) 

Intarfaoa  (•  publio  daolarations  •) 

Usas  owarlay 
,  globals 


prooadura  abort  ; 

procadura  add_toJC  (  tara  :  float  ;  st_row,  at_col  :  word  )  ; 
procadura  allooata_coluan  (  i  :  intagar  ;  r  :  word  )  ; 
prooadura  allooata_haight  ; 
prooadura  allooatajt  ; 

prooadura  allooatajvactor  (  VAX  X  :  ooluanjptr  )  ; 
prooadura  olaar_K  ; 

function  diagonal_araraga  (  VAX  darg  :  float  )  :  float  ; 
prooadura  diagonal_axtrsaas  (  VAX  daax,  dain  :  float  )  ; 
prooadura  KtiaasX  (  V,  X  :  ooluan_ptr  }  ; 
function  ain  (  a,  b  :  intagar  )  :  intagar  ; 
function  noraji  (  VAX  wnora  :  float  )  :  float  ; 
prooadura  print Jl  ; 

prooadura  ratriawa_X  {  VAX  h  :  aatrix_fila  )  ; 
prooadura  aolarr  (  n  :  intagar  )  ; 

prooadura  sat_r actor  (  v  :  ooluan  _ptr;  n  :  word;  v  :  float  )  ; 

procadura  stora_K  (  VAX  h  :  aatrix_fila  >  ; 

laplaaantation  (•  privata  daolarations  •) 


A-13 


(••••• 

•.  Abort  -  release  memory  And  di«. 
•) 

procedure  Abort  ; 
begin  («  Abort  •) 

raloAs# (  origin  )  ; 
holt (  l  ) 

•nd  (•  Abort  •)  ; 


•.  Add_to_K  -  Add  'tarn'  At  'at_row,at_col' 

•  where  at_row,  At_eol  give  tha  location  in 

•  A  "full"  matrix, 
a 

•) 

procadura  add_to_K  (  tara  :  float  ;  at_row.  at_col  :  word  )  ; 

VAR 

i  :  word  ; 

begin  (•  add_to_K  •) 

if  at_row  <■  at_col  than  begin  (•  K  ia  upper  triangle  only  •) 

i  :«  Suoo(  at_col  -  at_row  )  ; 

if  i  <«  height* [at_ool]  than  (•  inauranoa  •) 

K-[at_eoll*[i]  K*tat_ooiriil  ♦  tara 
and  (•  if  •) 

and  (•  add_to_K  •)  ; 


*.  alloc at a_coluan  -  arrange  storage  for  tha  i-th  column. 

a 

•) 

procadura  allocate_coluan  (  i  :  integer  ;  (•  tha  ooluan  of  interest  •) 

r  :  word  (•  the  sisa  desired  •) 

)  ; 

VAR 

o  :  longint  ; 
begin 

o  :■  r  •  sixeof.float  ; 

if  e  >  MeaArail  than  begin 

writaln  (  'Out  of  aeaory  attempting  to  allooata  o, 

'  bytes  to  ooluan  i,  '  -  allocata_coluan. '  )  ; 

abort 

and 

also  begin 

Oataaa  (  R*[il.  o  )  ; 
end 

end  (•  allocate_ooluan  *)  ; 
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•lloe*t«_h*ifht  -  allocate  storaf*  for  tho  column  haifhta. 


•> 

procedure  allocete_height  ; 

VAX 

e  :  lonfint  ; 

begin  (•  allocate.height  •) 

o  : »  neq  •  sizeof(  word  )  ; 

if  c  >  MeaAvail  then  begin 

writeln(  'Out  of  aeaory  etteepting  to  allocate  ',  c, 

*  byte*  to  eoluan  heights  -  allocate_height . '  ) 

abort 

end 

else 

Cetaea  (  height,  o  ) 
end  (•  allocate_height  *)  ; 


*.  allocate_K  -  allocate  storage  for  the  aatrix  of  coefficients. 
•) 

procedure  allocatejt  ; 

VAX 

c  :  longint  ; 
i  :  word  ; 

begin  (•  al locate Jt  •) 

c  :•  neq  *  sixeof(  pointer  )  ; 

if  c  >  MeaAvail  then  begin 

writeln  (  'Out  of  aeaory  atteapting  to  allocate  c, 

'  bytes  to  (K)  -  allocate.*. '  )  ; 

abort 

end 

else 

Cetaea  (  K,  o  )  ; 

for  i  :•  l  to  neq  do 

allocate.coluan (  i,  height* [i)  ) 

end  <•  allocatejt  •)  ; 
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allocat*_v*ctor  -  allocate  storage  for  a  'right-hand-* id*- like'  rector. 


I  T 

*. 

•) 

prooadur*  allooat*_v*ctor  {  VAX  X  :  coluan_ptr  )  ; 

VAX 

o  :  longlnt  ; 

begin  (*  alloeate_veetor  •) 

o  :■  neq  •  eizeof_float  ; 

if  c  >  HeaAvail  then  begin 

writeln  (  'Out  of  aeaory  attempting  to  allocate  c, 

'  byte*  to  coluan  rector  -  allooate_rector. '  )  ; 

abort 

end 

els* 

Getaea  (  X.  c  ) 
end  (•  allocate_r*ctor  •)  ; 


*.  cl*ar_K  -  act  the  contents  of  K  to  0.0. 
* 

•) 

procedure  cl*ar_K  ; 

VAX 

i  :  word  ; 
begin  (•  cl*ar_K  •) 

for  i  :«  l  to  neq  do 

s*t_rector(  K*[i],  height' JiJ,  0.0  ) 

end  (•  clear  Jt  *)  ; 


*.  diagonal_aeerage  -  return  arerage  of  ralues  on  diagonal  of  K. 


function  diagonal.arorage  (  VAX  darg  :  float  )  :  float  ; 
VAX 

i  :  word  ; 

begin  (•  diagonal_arerag*  •) 

darg  :*  0.0  ; 
for  i  :«  1  to  neq  do 
darg  darg  *  K'Ul'lU  S 
darg  :•  darg  /  neq  ; 
d i agona l_ar crag*  :•  darg 


md  (•  diagenaljarerag*  •)  ; 


•.  diagonal_sxtrsass  -  return  maximum  and  ainiaua  values  on  diagonal  of  R. 


procedure  diagonal_ext reams  (  VAX  daax,  da in  :  float  )  ; 
VAX  i  :  word  ; 
begin  (•  diagonal.extreaes  *) 
daax  K*[l]*[l]  ; 
dain  :•  R'tl]'|l]  ; 
for  i  :■  l  to  neq  do  begin 

if  R'tiJ'll]  >  daax  then  daax  R“[i]‘[l)  ; 
if  R‘li]'|l]  <  dain  then  dain  :■  R*[i)*[l]  ; 

end 

end  (•  diagonal_extreaes  •)  ; 


KtiaeaX  -  column  vector  F  :»  skyline  aatrix  X  tines  coluan  vector  X. 
uses  global  variables  'neq',  'K'  and  'height' 

K  is  syaaetrio: 
kll  k44 

k2l  k43 

k31  k42 
(sya)  k41 

Fart  1  uses  eaoh  coluan  froa  the  diagonal  to  the  aaxiaua 
height  in  eaoh  coluan  as  a  row.  This  ensures  that 
the  terns  to  the  left  of  the  diagonal  are  used  in  the 
aultiplioation. 

Part  2  adjusts  the  previous  suas  by  using  the  terns  above 
the  diagonal  as  the  reaainder  of  the  row. 


ocedure  RtiaesX  (  F,  X  :  coluan_ptr  )  ; 

VAX 

i.  j.  q  ••  word  ; 
s  :  float  ; 
begin  (*  RtiaesX  •) 

for  i  :■  1  to  neq  do  begin 
(•  Part  1  •) 
q  :■  height'll]  ; 
s  :<*  0.0  ; 

for  j  :■  1  to  height'll]  do  begin 
s  S  ♦  R‘[i]'[q]  •  X'UJ  : 
dec{  q  ) 
end  ; 

F'lU  s  ; 

(•  Fart  2  •) 

q  :»  i  -  1  ; 

for  j  :•  2  to  height'll]  do  begin 
F'tq]  F*|q]  ♦  K*|i]*|j]  •  **[1]  ! 
dec  (  q  ) 


wid  (•  for  •)  ; 
end  (•  RtiaesX  •)  ; 
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(••••* 

•  .  ain  -  ratumi  aaallaat  intagar  valua. 

•) 

function  ain  (  a,  b  :  intagar  )  :  intagar  ; 
begin  (•  ain  •) 

if  a  <■  b  than  ain  : -  a  alsa  ain  : -  b 
and  (•  ain  •)  ; 


e 

*.  norajt  -  a  nora  of  K. 

• 

•) 

function  norajt  (  VAR  wnora  :  float  )  :  float  ; 

VAR 

i  :  word  ; 
bagin  (•  norajt  •) 

▼nora  K*tl]*(l]  ; 
for  i  :■  1  to  neq  do 

if  K*[il*tU  >  ▼nora  than  ▼nora  K*UJ*U]  *. 
norajt  :»  vnora 
and  (•  norajt  •)  ; 


(•••*• 

* 

*.  print Jt  -  print  tha  ooaffieiant  aatrix  (for  dabugging  purpoaaa) 

•) 

procedure  print JC  ; 

VAR 

n  :  word  ; 
p  :  word  ; 
bagin  (•  printjt  •) 

for  n  :■  1  to  naq  do  bagin 

for  p  :■  l  to  haight*(n]  do  bagin 
writa<  K* [n]*(p]  ) 
and  (*  p  •)  : 
writaln 
and  (•  n  •) 


and  (•  printjt  •)  ; 
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retrieve_ K  -  fro*  <h>. 


• . 

*) 

procedure  r#tri#v#_K  (  VAR  h  :  aatrix_file  )  ; 
VAR 

i,  j  :  word  ; 
x  :  item  ; 

begin  (*  retrieve_K  •)' 

read (  h,  x.f  )  ; 
neq  :■  x.w  ; 

for  i  :■  1  to  neq  do  begin 

read(  h,  x.f  )  ; 
height* [i]  :»  x.w 
end  ; 

for  i  :■  1  to  neq  do 

for  j  l  to  height* ti]  do 
reed(  h.  K*[il*[j]  ) 

end  (•  retrievejk  •)  ; 


(•*•** 

•.  eet_vector  -  eet  V*[l..nJ  to  W. 

•) 

procedure  «et_vector  (  V  :  ooluxn_.pt r;  n  :  word;  V  :  float  )  ; 
begin  (•  eet_vector  •) 

for  n  :■  l  to  n  do 
V*tn]  :•  w 

end  (•  set_vector  *)  ; 


(*»•** 

•.  eolerr  -  print*  error  aesaage  and  die*. 

•) 

procedure  eolerr  (  n  :  integer  )  ; 
begin  (•  aolerr  •) 

writeln(  Diagonal  tera  in  ooluan  n,  '  ■  K*[n]*(l]  )  ; 

writeln(  '  Coefficient  aatrix  is  NOT  positive  definite.  -  solerr.'  ) 
halt(  l  ) 

end 


(•  solerr  •)  ; 


*.  «tora_K  -  im  K  on  <h>. 


procadura  stora_K  (  VAR  h  :  a*tri*_fila  )  ; 
VAR 

i,  j  :  word  ; 
x  :  itaa  ; 
basin  (•  atoraJC  *) 

x.w  :»  naq  ;  writa(  h,  x.t  )  ; 

for  i  :■  1  to  naq  do  basin 

x.w  :»  haisht*[i]  ; 
writa(  h.  x.f  ) 
and  ; 

for  i  :■  l  to  naq  do 

for  i  1  to  haight*[i]  do 

writa(  h.  K*[i]‘[j]  ) 

and  (*  atora_K  •)  ; 


and  (•  kxrutil  •) . 
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